program define summary_stats_forecasts, eclass

	syntax varlist(min=2 max=2 numeric), [K(int 4)] 
	
	local forecast_error: word 1 of `varlist'
	local actual: word 2 of `varlist'
		
	if "`method'" == "" local method "newey"
	
	tempvar MSE RMSE N mean_error mean_absolute_error MASE abs_error ///
		abs_error_naive mean_abs_error_naive res std_dev_x mean_x 
	
	* Create matrix for storing results
	matrix `res' = J(1, 10, .)
	matrix colnames `res' = N mean_x std_dev_x rho_x mean_error ///
				RMSE MASE rho_error adj_R_sq K
	
	* Calculate measures of forecast accuracy
	quietly {
		egen `MSE' = mean(`forecast_error' ^ 2)
		gen `RMSE' = sqrt(`MSE')
		egen `N' = count(`forecast_error')
		egen `mean_error' = mean(`forecast_error')
		egen `std_dev_x' = sd(`actual')
		egen `mean_x' = mean(`actual')
		gen `abs_error' = abs(`forecast_error')
		egen `mean_absolute_error' = mean(`abs_error')
		gen `abs_error_naive' = abs(`actual' - l.`actual')
		egen `mean_abs_error_naive' = mean(`abs_error_naive')
		gen `MASE' = `mean_absolute_error' / `mean_abs_error_naive'
	}

	* Calculate persistence of underlying process
	quietly reg `actual' l.`actual'
	
	* Collect results
	matrix `res'[1, 1] = `N'[1] 
	matrix `res'[1, 2] = `mean_x'[1]
	matrix `res'[1, 3] = `std_dev_x'[1]
	matrix `res'[1, 4] = _b[l.`actual']
	matrix `res'[1, 5] = `mean_error'[1]
	matrix `res'[1, 6] = `RMSE'[1]
	matrix `res'[1, 7] = `MASE'[1]
	
	* Calculate persistence of forecast errors
	quietly reg `forecast_error' l.`forecast_error'
	matrix `res'[1, 8] = _b[l.`forecast_error']
	
	* Estimate a regression of forecast errors on k
	* past forecast errors, get adjusted R-squared
	quietly {
	forvalues ll = 1/`k' {
		tempvar `forecast_error'_LAG`ll'
		gen ``forecast_error'_LAG`ll'' = l`ll'.`forecast_error'
	} 
	}
	
	local controls
	forvalues ll = 1/`k' {
		local controls `controls' ``forecast_error'_LAG`ll''
	}		
	
	quietly reg `forecast_error' `controls'
	matrix `res'[1, 9] = e(r2_a)
	matrix `res'[1, 10] = `k'

	ereturn matrix s `res'
		
end 
